/* ----------------------------------------------------------- */
/*                                                             */
/*                          ___                                */
/*                       |_| | |_/   SPEECH                    */
/*                       | | | | \   RECOGNITION               */
/*                       =========   SOFTWARE                  */
/*                                                             */
/*                                                             */
/* ----------------------------------------------------------- */
/* ----------------------------------------------------------- */
/*         Copyright:                                          */
/*                                                             */
/*         2011-2012  NUS                                      */
/*                    School of Computing                      */
/*                                                             */
/*   Use of this software is governed by a License Agreement   */
/*    ** See the file License for the Conditions of Use  **    */
/*    **     This banner notice must not be removed      **    */
/*                                                             */
/* ----------------------------------------------------------- */
/*        File: dspr.c: Discriminatave (MMI/MPE)  training */
/*     Using Frame Discrimination.                             */
/* ----------------------------------------------------------- */

#define EXITSTATUS 0 /*2 for gprof.*/

char *dspr_version = "!HVER!dspr:   3.4.1 [CUED 12/03/09]";
char *dspr_vc_id = "$Id: dspr.c,v 1.1.1.1 2006/10/11 09:55:01 jal58 Exp $";


#include "HShell.h"     /* HMM ToolKit Modules */
#include "HMem.h"
#include "HMath.h"
#include "HSigP.h"
#include "HAudio.h"
#include "HWave.h"
#include "HVQ.h"
#include "HParm.h"
#include "HLabel.h"
#include "HModel.h"
#include "HTrain.h"
#include "HUtil.h"
#include "HAdapt.h"
#include "HFB.h"
#include "HDict.h"
#include "HNet.h"
#include "HLM.h"
#include "HLat.h"
#include "HArc.h"
#include "HFBLat.h"
#include "HExactMPE.h"

/* Parameters Global storage  */


#define MAX(a,b) ((a)>(b)?(a):(b))
#define MIN(a,b) ((a)<(b)?(a):(b))
#define ABS(a) ((a)>0?(a):-(a))
#define FINITE(x) (!isnan(x) && x<1.0e+30 && x>-1.0e+30)


/* Trace Flags */
#define T_TOP   0001    /* Top level tracing */
#define T_TIM   0002    /* Output timings */

/* possible values of updateMode */
#define UPMODE_DUMP 1
#define UPMODE_UPDATE 2
#define UPMODE_BOTH 3


int nDenLats = 0;
int nNumLats = 0;

static char *latFileMask = NULL;
static char *LatMask_Numerator   = NULL;
static char *LatMask_Denominator = NULL;

/* Global Settings */

static char * hmmDir = NULL;     /* directory to look for hmm def files */
static char * hmmExt = NULL;     /* hmm def file extension */
static char * newDir = NULL;     /* directory to store new hmm def files */
static char * newExt = NULL;     /* extension of new reestimated hmm files */

static char * statFN=0;            /* stats file, if any */
static float minVar  = 0.0;      /* minimum variance (diagonal only).
                                    This gives a constant floor to all dimensions' vars. */

static int updateMode = UPMODE_UPDATE; /* dump summed accs, update models or do both? */

static float mixWeightFloor=MINMIX*2; /* Floor for mixture weights */
/* There *was* an initialisation of 0.0.  However, this interacts very badly
   with the discriminative technique, which really relies on them being
   taken back up (? I wrote this earlier, dont know if it is still true)
   I have changed it further from 1.01 to 2, which I had usually put in
   the HTE file.  This enables less stuff in the HTE without me having to
   do experiments to determine the best value.
*/

static UPDSet uFlags = UPMEANS|UPVARS|UPTRANS|UPMIXES;   /* update flags */
static UPDSet uFlagsAccs = UPMEANS|UPVARS|UPTRANS|UPMIXES;   /* used in storing accs. */
static UPDSet uFlagsMLE = 0; /*which we only update with MLE, ignoring the MMI parameters.*/


static int parMode   = 0;       /* enable one of the parallel modes */
/* i.e.  0 for reestimation, 1,2,3... for accumulating .acc files. */


static Boolean stats = FALSE;    /* enable statistics reports */
static char * mmfFn  = NULL;     /* output MMF file, if any */
static long int trace     = 1;        /* Trace level */

static Boolean saveBinary = FALSE;  /* save output HMMs in binary  */

static Boolean ML_MODE = FALSE;     /* when only one set of accs are supplied. */
static Boolean MPE = FALSE;         /* when we are doing MPE/MWE. */
static Boolean MMIPrior = FALSE;    /* use MMI prior as I-smoothing prior */
static float MMITauI = 0.0;        /* I-smoothing tau for MMI prior in MPE training */

static float E = 2.0;               /* constant used in BW updatel */

static float CWeights = 1.0;
static float CTrans = 1.0;


static Boolean MPEStoreML=FALSE;   /*  Set TRUE if we need to accumulate ML stats while doing MPE. */
static Boolean THREEACCS = FALSE;  /*  Set TRUE if 3 sets of accs need to be stored (for MPE, or possibly for MMI-MAP) */
static int NumAccs;                /*  Set in Initialise() to  1 or 2 or 3. */

static float ISmoothTau = 0.0;        /* I-smoothing: a h-crit-like thing.  Set to 100 for MMI or 50 for MPE, or 25 for MWE. */

/* If unset and ISmoothTau is set, try setting these to 1/10 the occupancy for ISmoothTau. */
static float ISmoothTauTrans = 0.0;
static float ISmoothTauWeights = 0.0;
static Boolean ISmoothTauTransSet=FALSE;
static Boolean ISmoothTauWeightsSet=FALSE;

/* R.E. priors: */
static float PriorTau = 0.0;         /* tau value [e.g. 10,25] for use in discriminative MAP with -Hprior option. */
static float PriorK = 0.0;           /* e.g. use 1 if disc. trained prior model is better than ML trained new model.  */

static float PriorTauWeights = 0.0;
static float PriorTauTrans = 0.0;


static Boolean twoDataFiles = FALSE; /* For training using two data files, probably never needed & code not tested. */

static float MinOcc = 10;              /* Minimum numerator (ML) occupancy for a Gaussian to be updated */
static float MinOccTrans = 10;         /* Minimum numerator (ML) occupancy for a transition row */
static float MinOccWeights = 10;       /* Minimum numerator (ML) occupancy for a set of weights. */

static float hcrit = 1.0; /* Scale on denominator (MMI) part. Not really useful.  */

static float varFloorPercent = 0;
static float varSmooth = 0;

static Boolean useLLF = FALSE;          /* use directory based LLF files instead of individual lattices */
/* Global non-config variables */

static ConfParam *cParm[MAXGLOBS];   /* configuration parameters */
static int nParm = 0;                /* total num params */

static int maxM = 0;             /* max mixtures in any model */
static int vSize;                /* input vector size */
static int S;                    /* number of data streams */   /*! Equals 1 or error![?] */
static HSetKind hsKind;          /* kind of the HMM system */ /*!Must be PLAINHS | SHAREDHS (|TIEDHS?) */

static int L;                        /* number of logical HMM's */
static int P;                        /* number of physical HMM's */
static HMMSet hset;                  /* Set of HMMs to be re-estimated */
static HMMSet hset2; /* store info of MMF2 (after SPR)*/
static HMMSet hset_prior;            /* Usually uninitialised, except for MMI-MAP/MPE-MAP */
Boolean hset_prior_initialised = FALSE;

static FBLatInfo  fbInfo;            /* Structure for discriminative forward-backward. */



static int totalT=0;                             /* total number of frames in training data */
static LogDouble totalPr1=0,totalPr2=0;              /* total log prob upto current utterance, totalPr3 for MMI den in MPE with MMI Prior*/


static Vector vFloor[SMAX];          /* variance floor - default is all zero */

static MemHeap accStack;           /* accumulated statistics */
static MemHeap hmmStack;           /* HMM defs and related structures */
static MemHeap transStack;         /* Transcriptions... comment in original HERest says transformations, but it's wrong.*/

static MemHeap latStack;           /* Lattices. */


static MemHeap hmmStack2;

Vocab vocab;

/* MPE: */
static int TotalNWords=0;
static double TotalCorr=0;

/* information about transforms */
static XFInfo xfInfo;

/* static prior */
static Boolean STATICPRIOR = FALSE;

/* Thang v110701 */
static char * accsPrefix = NULL;     /* to look for files accsPrefix(1|2|3).(info|end|mean|mean.occ|var|var.occ|trace|trace.occ|weight|weight.occ) */

/* possible values of accsMode */
#define ACCS_ASCII_IN 1
#define ACCS_ASCII_OUT 2
static int accsMode = 0;     /* when only one set of accs are supplied. */

#define T_ALD  00400       /* Trace acc load/dump */
static Boolean inLdBinary = FALSE;       /* load in ascii/binary */
static Boolean outLdBinary = FALSE;       /* dump in ascii/binary */

void baoUpdateModels(DVector * mean,DVector *var,int numFeatures);


/* ------------------ Process Command Line -------------------------- */

/* SetConfParms: set conf parms relevant to dspr  */

void SetConfParms(void)
{
   int i;
   Boolean b;
   double f;
   char buf[MAXSTRLEN];

   nParm = GetConfig("dspr", TRUE, cParm, MAXGLOBS);
   if (nParm>0) {
      if (GetConfInt(cParm,nParm,"TRACE",&i)) trace = i;
      if (GetConfFlt(cParm,nParm,"VARFLOORPERCENTILE",&f)) varFloorPercent = f;
      if (GetConfFlt(cParm,nParm,"VARSMOOTH",&f)) varSmooth = f;
      if (GetConfFlt(cParm,nParm,"C",&f)){ CWeights = CTrans = f; }
      if (GetConfFlt(cParm,nParm,"CW",&f)) CWeights = f;
      if (GetConfFlt(cParm,nParm,"CT",&f)) CTrans = f;
      if (GetConfFlt(cParm,nParm,"MINOCC",&f)){ MinOcc = f; }
      if (GetConfFlt(cParm,nParm,"MINOCCTRANS",&f)){ MinOccTrans = f; }
      if (GetConfFlt(cParm,nParm,"MINOCCWEIGHTS",&f)){ MinOccWeights = f; }
      if (GetConfBool(cParm,nParm,"SAVEBINARY",&b)) saveBinary = b;
      if (GetConfFlt(cParm,nParm,"E",&f)) E = f;
      if (GetConfFlt(cParm,nParm,"DFACTOROCC",&f)) E = f; /*Back-compat. */
      if (GetConfFlt(cParm,nParm,"HCRIT",&f)) hcrit = f;
      if (GetConfBool(cParm,nParm,"MPE",&b)){  MPE = b; THREEACCS=MPE&&MPEStoreML; }
      if (GetConfBool(cParm,nParm,"MWE",&b)){  MPE = b; THREEACCS=MPE&&MPEStoreML; } /* "MWE" has identical effects here but differs in HFBLat.c */
      if (GetConfBool(cParm,nParm,"MEE",&b)){  MPE = b; THREEACCS=MPE&&MPEStoreML; } /* Back-compat. */
      if (GetConfBool(cParm,nParm,"MLE",&b)){  ML_MODE=TRUE; THREEACCS=FALSE;
      uFlagsMLE =  UPMEANS|UPVARS|UPTRANS|UPMIXES; }
      if (GetConfBool(cParm,nParm,"MMIPRIOR",&b)){  MMIPrior = b;}
      if (GetConfFlt(cParm,nParm,"MMITAUI",&f)){ MMITauI = f;}
      if (GetConfFlt(cParm,nParm,"ISMOOTHTAU",&f)){ ISmoothTau = f; MPEStoreML=TRUE; THREEACCS=MPE&&MPEStoreML; }
      if (GetConfFlt(cParm,nParm,"ICRITOCC",  &f)){ ISmoothTau = f; MPEStoreML=TRUE; THREEACCS=MPE&&MPEStoreML; } /*back-compat. */

      if (GetConfFlt(cParm,nParm,"ISMOOTHTAUT",&f)){ ISmoothTauTrans = f; ISmoothTauTransSet=TRUE; }
      if (GetConfFlt(cParm,nParm,"ISMOOTHTAUW",&f)){ ISmoothTauWeights = f; ISmoothTauWeightsSet=TRUE; }

      if (GetConfFlt(cParm,nParm,"PRIORTAU",&f)){ PriorTau = f; }
      if (GetConfFlt(cParm,nParm,"PRIORTAUW",&f)){ PriorTauWeights = f; }
      if (GetConfFlt(cParm,nParm,"PRIORTAUT",&f)){ PriorTauTrans = f; }
      if (GetConfFlt(cParm,nParm,"PRIORK",&f)){ PriorK = f; }
      if (GetConfBool(cParm,nParm,"STATICPRIOR",&b)) {
         STATICPRIOR=b;
         if (STATICPRIOR) PriorK = 1.0; else PriorK = 0.0;
      }

      if (GetConfFlt(cParm,nParm,"MIXWEIGHTFLOOR",&f)){ mixWeightFloor = MINMIX * f; }


      if (GetConfStr(cParm,nParm,"LATFILEMASK",buf)) {
         latFileMask = (char*)malloc(strlen(buf)+1);
         strcpy(latFileMask, buf);
      }

      if (GetConfStr(cParm,nParm,"LATMASKNUM",buf)) {
         LatMask_Numerator = (char*)malloc(strlen(buf)+1);
         strcpy(LatMask_Numerator,buf);
      }
      if (GetConfStr(cParm,nParm,"LATMASKDEN",buf)) {
         LatMask_Denominator = (char*)malloc(strlen(buf)+1);
         strcpy(LatMask_Denominator,buf);
      }
      if (GetConfStr(cParm,nParm,"INXFORMMASK",buf)) {
         xfInfo.inSpkrPat = (char*)malloc(strlen(buf)+1);
         strcpy(xfInfo.inSpkrPat,buf);
      }
      if (GetConfStr(cParm,nParm,"PAXFORMMASK",buf)) {
         xfInfo.paSpkrPat = (char*)malloc(strlen(buf)+1);
         strcpy(xfInfo.paSpkrPat,buf);
      }

      if (GetConfBool(cParm,nParm,"USELLF",&b))  useLLF = b;

      if (GetConfStr(cParm,nParm,"UPDATEMODE",buf)) {
         if (!strcmp (buf, "DUMP")) updateMode = UPMODE_DUMP;
         else if (!strcmp (buf, "UPDATE")) updateMode = UPMODE_UPDATE;
         else if (!strcmp (buf, "BOTH")) updateMode = UPMODE_BOTH;
         else HError(2319, "Unknown UPDATEMODE specified (must be DUMP, UPDATE or BOTH)");
      }

			/* Thang v110701: binary option */
			if (GetConfBool(cParm,nParm,"INBINARYACCFORMAT",&b)) inLdBinary = b;
			if (GetConfBool(cParm,nParm,"OUTBINARYACCFORMAT",&b)) outLdBinary = b;

   }

   if (MPE && uFlagsMLE) HError(1, "Can't combine MPE with ML update of some parameters (code could be simply added).");
   if (MMIPrior && !THREEACCS) HError(999, "MMI Prior must be used in MPE update (THREEACCS).");
   if ((STATICPRIOR && PriorK==0.0) || (!STATICPRIOR && PriorK==1.0)) HError(999, "Specify either PRIORK or STATICPRIOR (PRIORK overwrites value given by STATICPRIOR).");
   /*   if(ISmoothTau && !ISmoothTauTransSet){ ISmoothTauTrans = 10; printf("Smoothing transitions with tau=%f since ISMOOTHTAUT not set\n",ISmoothTauTrans); }
        if(ISmoothTau && !ISmoothTauWeightsSet){ ISmoothTauTrans = 10; printf("Smoothing weights with tau=%f since ISMOOTHTAUW not set\n",ISmoothTauWeights); } */
}


void ReportUsage(void)
{
   printf("\nUSAGE: dspr [options] hmmList dataFiles...\n\n");
   printf(" Option                                   Default\n\n");
	 /* Thang v110701: add options */
   printf(" -accs s     Prefix to name input/output accs files\n");
   printf(" -asciiin    Expected additional accs files in ascii (compensated mean & variance)\n");
   printf(" -asciiout   Extract mean and variance to be outputed to ascii file");
   printf(" -update s   Update mode (DUMP, UPDATE, BOTH).                  UPDATE   \n");
   printf(" -inB    Input accs files (HDR*) are expected to be in binary format                 off\n");
   printf(" -outB   Output accs files (HDR*) will be dumped in binary format                 off\n");
   printf(" -H  s path to MMF file 1.\n");
   printf(" -P  s path to MMF file 2.\n");
   printf(" -M  s path to save the ouput file.\n");
   printf(" -C  s Configure file.\n");
   printf(" -p N    set parallel mode to N             off\n");
   PrintStdOpts("BEFHIMSTJX"); /*K,G,L removed*/
   printf(" Note: doesn't work if means and variances are shared independently.\n");
   printf("\n\n");
}

void SetuFlags(UPDSet *uFlags)
{
   char *s;

   s=GetStrArg();
   *uFlags=0;
   while (*s != '\0')
      switch (*s++) {
      case 't': (*uFlags)+=UPTRANS; break;
      case 'm': (*uFlags)+=UPMEANS; break;
      case 'v': (*uFlags)+=UPVARS;  break;
      case 'w': (*uFlags)+=UPMIXES; break;
      case 'a': (*uFlags)+=UPXFORM; break;
      }
}


void PrintCriteria(){
   printf("\nMMI criterion per frame is: %f (%f - %f)\n", (totalPr1-totalPr2)/totalT,totalPr1/totalT,totalPr2/totalT);
   if(MPE) printf("\nMPE/MWE criterion is: %f ( %f / %d )\n", TotalCorr/TotalNWords, TotalCorr, TotalNWords);
   if(!MPE || MPEStoreML) printf("\nML criterion per frame is: %f (%f/%d)\n", totalPr1/totalT,totalPr1, totalT);
}

/*******************************************/
/********* CONVERT BINARY TO ASCII *********/
/*******************************************/
/* LoadWtAcc: new inc of wt acc from file f */
static void LoadWtAcc(Source *src, WtAcc *wa, int numMixtures, Boolean ldBinary)
{
   int m;
   float f;
   Vector cTemp;

   cTemp = CreateVector(&gstack,numMixtures);
   ReadVector(src,cTemp,ldBinary);
   for (m=1;m<=numMixtures;m++){
     if(!finite(cTemp[m]))
       HError(7191, "Infinite WtAcc!");
     wa->c[m] += cTemp[m];
   }
   ReadFloat(src,&f,1,ldBinary);
   wa->occ += f;
   FreeVector(&gstack,cTemp);
}

/* LoadMuAcc: new inc of mean acc from file f */
static void LoadMuAcc(Source *src, MuAcc *ma, int vSize, Boolean ldBinary)
{
   int k;
   Vector vTemp;
   float f;

   vTemp = CreateVector(&gstack,vSize);
   ReadVector(src,vTemp,ldBinary);
   for (k=1;k<=vSize;k++){
     if(!finite(vTemp[k]))
       HError(7191, "Infinite MuAcc!");
      ma->mu[k] += vTemp[k];
   }
   ReadFloat(src,&f,1,ldBinary);
   ma->occ += f;
   FreeVector(&gstack,vTemp);
}

/* LoadVaAcc: p'th inc of variance acc from file f */
static void LoadVaAcc(Source *src, VaAcc *va, int vSize, CovKind ck, Boolean ldBinary)
{
   int k,kk;
   Vector vTemp;
   TriMat mTemp;
   float f;

   switch(ck){
   case DIAGC:
   case INVDIAGC:
      vTemp = CreateVector(&gstack, vSize);
      ReadVector(src,vTemp,ldBinary);
      for (k=1;k<=vSize;k++){
        if(!finite(vTemp[k]))
           HError(7191, "Infinite VaAcc!");
         va->cov.var[k] += vTemp[k];
      }
      FreeVector(&gstack, vTemp);
      break;
   case FULLC:
   case LLTC:
      mTemp = CreateTriMat(&gstack,vSize);
      ReadTriMat(src,mTemp,ldBinary);
      for (k=1;k<=vSize;k++)
         for (kk=1; kk<=k; kk++) {
            va->cov.inv[k][kk] += mTemp[k][kk];
         }
      FreeTriMat(&gstack,mTemp);
      break;
   }
   ReadFloat(src,&f,1,ldBinary);
   va->occ += f;
}

/* LoadTrAcc: p'th inc of transition acc from file f */
static void LoadTrAcc(Source *src, TrAcc *ta, int numStates, Boolean ldBinary)
{
   int i,j;
   Matrix tTemp;
   Vector nTemp;

   tTemp = CreateMatrix(&gstack,numStates,numStates);
   nTemp = CreateVector(&gstack,numStates);
   ReadMatrix(src,tTemp,ldBinary);
   for (i=1;i<=numStates;i++)
      for (j=1;j<=numStates;j++)
         ta->tran[i][j] += tTemp[i][j];
   ReadVector(src,nTemp,ldBinary);
   for (i=1;i<=numStates;i++)
      ta->occ[i] += nTemp[i];
   FreeMatrix(&gstack,tTemp);
}

/* CheckPName: check dumped name matches hmm phys name */
static void CheckPName(Source *src, char *pname)
{
   int c;
   char buf[MAXSTRLEN];

   ReadString(src,buf);
   c = GetCh(src);
   if (c != '\n')
      HError(7150,"CheckPName: Cant find EOL");
   if (strcmp(pname,buf) != 0)
      HError(7150,"CheckPName: expected %s got %s",pname,buf);
}

/* CheckMarker: check file f has a marker next */
static void CheckMarker(Source *src,Boolean ldBinary)
{
   int mark = 123456, temp;

   ReadInt(src,&temp,1,ldBinary);
   if (temp != mark)
      HError(7150,"CheckMarker: Marker Expected in Dump File");
}

/* GetDumpFile: Process dump file name and open it */
static FILE * GetDumpFile(char *name, int n)
{
   char buf[MAXSTRLEN],num[20];
   int i,j,k,len,nlen;
   FILE *f;

   sprintf(num,"%d",n);
   len = strlen(name); nlen = strlen(num);
   for (i=0,j=0; i<len; i++) {
      if (name[i] == '$')
         for (k=0; k<nlen; k++)
            buf[j++] = num[k];
      else
         buf[j++] = name[i];
   }
   buf[j] = '\0';
   f = fopen(buf,"wb"); /* Binary file */
   if (f==NULL)
      HError(7111,"GetDumpFile: cannot open acc dump file %s",buf);
   if (trace & T_ALD)
      printf("Dumping accumulators to file %s\n",buf);
   return f;
}

static void WriteVectorToAscii(FILE * f, DVector v, int numElements){
    int i = 1;
    for (; i< numElements; i++){
      fprintf(f, "%e ", v[i]);
    }
    fprintf(f, "%e\n", v[numElements]);
}

static void printVector(DVector v, int numElements){
    int i = 1;
    for (; i<= numElements; i++){
      fprintf(stderr, "%e ", v[i]);
    }
    fprintf(stderr, "\n");
}

static float sumVector(DVector v, int numElements){
    int i = 1;
		float sum = 0;
    for (; i<= numElements; i++){
			sum += abs(v[i]);
    }
		return sum;
}

static void printMatrix(Matrix m, int numElements){
	int i = 1;
	for (; i<= numElements; i++){
		int j = 1;
		for (; j<= numElements; j++){
			fprintf(stderr, "%e ", m[i][j]);
		}
		fprintf(stderr, "\n");
	}
}

/*LoadMMF: get mean and var in hset*/
void loadMMF(HMMSet *hset, Boolean ldBinary,DVector *stateNoiseMeanData,DVector *stateNoiseVarData, char **arNoiseInfo)
{

	HLink hmm;
	HMMScanState hss;
	int size;
	MixPDF* mp;

	int numComponents = 0; int numStates = 0; 	int numHMMs = 0;

	int numFeatures = hset->vecSize;

	int numMix = hset->numMix;

	DVector vTemp = NULL;
	vTemp = CreateDVector(&gstack,numFeatures);
	int ii = 1;
	NewHMMScan(hset, &hss);
	do {
		numHMMs++;
		hmm = hss.hmm;

		int stateId = 1;
		while (GoNextState(&hss,TRUE)) {
			stateId++;

			while (GoNextStream(&hss,TRUE)) {
				size = hset->vecSize;

				if (hss.isCont){

					int numMixtures = 0;
					numStates++;
					while (GoNextMix(&hss,TRUE)) {
						numMixtures++;
						numComponents++;

						TouchV(hss.mp->mean);

						for(ii=1; ii<=numFeatures; ii++)
							vTemp[ii] = hss.mp->mean[ii];

						stateNoiseMeanData[numComponents] = CreateDVector(&gstack,numFeatures);
						CopyDVector(vTemp,stateNoiseMeanData[numComponents]);
						TouchV(hss.mp->cov.var);

						for(ii=1; ii<=numFeatures; ii++){

							vTemp[ii] = hss.mp->cov.var[ii];
						}

						stateNoiseVarData[numComponents] = CreateDVector(&gstack,numFeatures);
						CopyDVector(vTemp,stateNoiseVarData[numComponents]);

					} /* while (GoNextMix(&hss,TRUE)) */

						char* phone = malloc(15*sizeof(char));
						char *strMixture = malloc(5*sizeof(char));;
						sprintf(phone, "h-%s-s%d", hss.mac->id->name, stateId);
						sprintf(strMixture, "%d", numMixtures);
						arNoiseInfo[3*numStates] = phone;
						arNoiseInfo[3*numStates+1] = strMixture;
						arNoiseInfo[3*numStates+2] = strMixture;

				} /* if (hss.isCont) */
			} /* while (GoNextStream(&hss,TRUE)) */
		} /*  while (GoNextState(&hss,TRUE)) */


	} while (GoNextHMM(&hss));
	EndHMMScan(&hss);
	/* load mean and var from tied*/

}

/* EXPORT->LoadAccs: inc accumulators in hset by vals in fname */
Source thangLoadAccsParallel(HMMSet *hset, char *fname, UPDSet uFlags, int index, char *outName,Boolean ldBinary,DVector *stateCleanMeanData,DVector *stateCleanVarData, char **arCleanInfo)
{
	Source src;
	HLink hmm;
	HMMScanState hss;
	int size,negs,m,s,phoneId=0;
	MixPDF* mp;

	int numComponents = 0; int numStates = 0; 	int numHMMs = 0;
	char inOutName[256]; sprintf(inOutName, "%s.%s", outName, "info"); /* HMM state, mixture + trace info */
	char muOutName[256]; sprintf(muOutName, "%s.%s", outName, "mean");
	char vaOutName[256]; sprintf(vaOutName, "%s.%s", outName, "var");
	FILE * inof = NULL; FILE * muof = NULL; FILE * vaof = NULL;

	int numFeatures = hset->vecSize;

	int numMix = hset->numMix;
	DVector vTemp = NULL;
	if(accsMode&ACCS_ASCII_OUT){
		vTemp = CreateDVector(&gstack,vSize);
		inof = GetDumpFile(inOutName, 0);
		muof = GetDumpFile(muOutName, 0);
		vaof = GetDumpFile(vaOutName, 0);
		fprintf(inof, "%010d\n", 0); /* placeholder for num components */
		fprintf(inof, "%d\n", numFeatures); /* feature dimension */
	}

	if (trace & T_ALD)
		printf("Loading accumulators from file %s\n",fname);

	if(InitSource(fname,&src,NoFilter)<SUCCESS)
		HError(7110,"LoadAccs: Can't open file %s", fname);
	NewHMMScan(hset, &hss);
	do {
		numHMMs++;

		hmm = hss.hmm;
		CheckPName(&src,hss.mac->id->name);

		ReadInt(&src,&negs,1,ldBinary);
		negs += (int)hmm->hook; hmm->hook = (void *)negs;

		int stateId = 1;
		while (GoNextState(&hss,TRUE)) {
			stateId++;

			while (GoNextStream(&hss,TRUE)) {
				if ((uFlags&UPSEMIT) && (strmProj)) size = hset->vecSize;
				else size = hset->swidth[hss.s];

				LoadWtAcc(&src,((WtAcc *)hss.ste->hook)+index,hss.M, ldBinary);

				if (hss.isCont){
					int numMixtures = 0;
					numStates++;

					while (GoNextMix(&hss,TRUE)) {
						numMixtures++;

						numComponents++;
						MuAcc* tmpMu = NULL;
						float muOcc = 0;
						if ((uFlags&UPMEANS) && (!IsSeenV(hss.mp->mean))) {

							LoadMuAcc(&src,((MuAcc *)GetHook(hss.mp->mean))+index,size, ldBinary);
							TouchV(hss.mp->mean);

							/* Thang v110701: mu */
							if(accsMode&ACCS_ASCII_OUT){
								tmpMu = ((MuAcc *)GetHook(hss.mp->mean))+index;

								/* compute accs mean */
								int ii = 1;
								muOcc = tmpMu->occ;
								if (muOcc < 1e-40) {
									for(; ii<=numFeatures; ii++){
										vTemp[ii] = 0;
									}
									/* WriteVectorToAscii(muof,tmpMu->mu,numFeatures); */
								} else {

									for(; ii<=numFeatures; ii++){

										vTemp[ii] = hss.mp->mean[ii] + tmpMu->mu[ii]/muOcc; /* mu_old + mu/occ */


									}

								}

								stateCleanMeanData[numComponents] = CreateDVector(&gstack,numFeatures);
								CopyDVector(vTemp,stateCleanMeanData[numComponents]);

								WriteVectorToAscii(muof,vTemp,numFeatures);
							}
						}

						if ((uFlags&UPSEMIT) && (!IsSeenV(hss.mp->cov.var))) {
							LoadVaAcc(&src,((VaAcc *)GetHook(hss.mp->cov.var))+index,
									size,FULLC,ldBinary);
							TouchV(hss.mp->cov.var);
						} else if ((uFlags&UPVARS) && (!IsSeenV(hss.mp->cov.var))) {

							LoadVaAcc(&src,((VaAcc *)GetHook(hss.mp->cov.var))+index,size,hss.mp->ckind,ldBinary);
							TouchV(hss.mp->cov.var);

							/* Thang v110701: va */
							if(accsMode&ACCS_ASCII_OUT){
								if (!(uFlags&UPMEANS)) {
									HError(-1, "dspr: can't compute new var in this condition since mu is undefined");
								}
								VaAcc* tmpVa = ((VaAcc *)GetHook(hss.mp->cov.var))+index;

								float vaOcc = tmpVa->occ;
								int ii = 1;
								if (muOcc < 1e-40 || vaOcc < 1e-40) {
									if(vaOcc > 0){
										/*fprintf(stderr, "# Small muOcc=%e vaOcc=%e, component %d\n", muOcc, vaOcc, numComponents);*/
									}

									for(; ii<=numFeatures; ii++){
										vTemp[ii] = 0;
									}
									/* WriteVectorToAscii(vaof,tmpVa->cov.var,numFeatures); */
								} else {
									/* compute accs var */
									for(; ii<=numFeatures; ii++){
										float ratio = tmpMu->mu[ii]/muOcc;
										vTemp[ii] = tmpVa->cov.var[ii]/vaOcc - ratio*ratio; /* var/occVar - (mu/occMu)^2 */
										if(fabs(vTemp[ii]) < 1e-3 && vTemp[ii] <= 0){
											if(fabs(vTemp[ii]) > 1e-4){
												/*fprintf(stderr, "# Medium negative variance=%e, component %d, id=%d, var/occ = %e/%e = % e - (mu/occ=%e/%e)^2=%e\n", vTemp[ii], numComponents, ii, tmpVa->cov.var[ii], vaOcc, tmpVa->cov.var[ii]/vaOcc, tmpMu->mu[ii], muOcc, ratio*ratio);*/
											}
											vTemp[ii] = 1e-10;
										} else if(vTemp[ii] < 0){
											/*fprintf(stderr, "# Large negative variance=%e, component %d, id=%d, var/occ = %e/%e = % e - (mu/occ=%e/%e)^2=%e\n", vTemp[ii], numComponents, ii, tmpVa->cov.var[ii], vaOcc, tmpVa->cov.var[ii]/vaOcc, tmpMu->mu[ii], muOcc, ratio*ratio);*/
											exit(1);
										}
									}

								} /* end if occ */
								stateCleanVarData[numComponents] = CreateDVector(&gstack,numFeatures);
								CopyDVector(vTemp,stateCleanVarData[numComponents]);
								WriteVectorToAscii(vaof,vTemp,numFeatures);
							} /* if(accsMode&ACCS_ASCII_OUT) */
						}
					} /* while (GoNextMix(&hss,TRUE)) */

					if(accsMode&ACCS_ASCII_OUT){
						fprintf(inof, "h-%s-s%d\n%d\n%d\n", hss.mac->id->name, stateId, numMixtures, numMixtures);
						char* phone = malloc(15*sizeof(char));
						char *strMixture = malloc(5*sizeof(char));;
						sprintf(phone, "h-%s-s%d", hss.mac->id->name, stateId);
						sprintf(strMixture, "%d", numMixtures);
						arCleanInfo[3*numStates] = phone;
						arCleanInfo[3*numStates+1] = strMixture;
						arCleanInfo[3*numStates+2] = strMixture;

					}
				} /* if (hss.isCont) */
			} /* while (GoNextStream(&hss,TRUE)) */
		} /*  while (GoNextState(&hss,TRUE)) */


		if (!IsSeenV(hmm->transP)){
			LoadTrAcc(&src, ((TrAcc *) GetHook(hmm->transP))+index,hss.N,ldBinary);
			TouchV(hmm->transP);
		}

		CheckMarker(&src,ldBinary);
	} while (GoNextHMM(&hss));
	EndHMMScan(&hss);

	/*fprintf(stderr, "Num states = %d, Num HMMs = %d, Num HMM components = %d\n", numStates, numHMMs, numComponents);*/

	int numTiedComponents = 0;
	if (hset->hsKind == TIEDHS){
		if(hset->swidth[0] > 1){
			HError(-1, "dspr: not handle multiple streams yet\n");
		}
		numTiedComponents = hset->tmRecs[0].nMix;
		numComponents += numTiedComponents;

		for (s=1; s<=hset->swidth[0]; s++){
			size = hset->swidth[s];
			for (m=1;m<=hset->tmRecs[s].nMix; m++){
				mp = hset->tmRecs[s].mixes[m];
				LoadMuAcc(&src,((MuAcc *)GetHook(mp->mean))+index,size,ldBinary);
				fprintf(stderr, "# Tied states\n");

				LoadVaAcc(&src,((VaAcc *)GetHook(mp->cov.var))+index,size,mp->ckind,ldBinary);

				if(accsMode&ACCS_ASCII_OUT){


					/* TODO: fprintf(inof, "h-%s-s%d\n%d\n%d\n", hss.mac->id->name, stateId, numMixtures, numMixtures); */

					MuAcc* tmpMu = ((MuAcc *)GetHook(hss.mp->mean))+index;
   				WriteVectorToAscii(muof,tmpMu->mu,numFeatures);

					VaAcc *	tmpVa = ((VaAcc *)GetHook(hss.mp->cov.var))+index;
      		WriteVectorToAscii(vaof,tmpVa->cov.var,numFeatures);
				}
			}
		}
	}

	if(accsMode&ACCS_ASCII_OUT){
	/*	fprintf(stderr, "Num tied components = %d\n", numTiedComponents);*/
		/*fprintf(stderr, "Num components = %d\n", numComponents);*/

		/* Add number of components at the beginning of file */
		fseek(inof, 0, SEEK_SET);
		fprintf(inof, "%d\n", numStates);

		/*FreeVector(&gstack,vTemp);*/

		/* close files */
		fclose(inof);
		fclose(muof);
		fclose(vaof);
	}
	return src;
}

/*******************************************/
/********* CONVERT ASCII TO BINARY *********/
/*******************************************/
/* DumpPName: dump physical HMM name */
static void DumpPName(FILE *f, char *pname)
{
   WriteString(f,pname,DBL_QUOTE);
   fprintf(f,"\n");
}

/* DumpMarker: dump a marker into file f */
static void DumpMarker(FILE *f, Boolean ldBinary) /* Thang v110701: add ldBinary */
{
   int mark = 123456;

   WriteInt(f,&mark,1,ldBinary);
   if (!ldBinary) fprintf(f,"\n");
}

/* DumpWtAcc: dump wt acc to file f */
static void DumpWtAcc(FILE *f, WtAcc *wa, Boolean ldBinary)
{
   WriteVector(f,wa->c,ldBinary);
   WriteFloat(f,&(wa->occ),1,ldBinary);
   if (!ldBinary) fprintf(f,"\n");
}

/* DumpMuAcc: dump mean acc to file f */
static void DumpMuAccFromAscii(FILE *f, Vector mu, float occ, Boolean ldBinary)
{
   WriteVector(f,mu,ldBinary);
   WriteFloat(f,&occ,1,ldBinary);
   if (!ldBinary) fprintf(f,"\n");
}

/* DumpVaAcc: dump variance acc to file f */
static void DumpVaAccFromAscii(FILE *f, Vector va, float occ, Boolean ldBinary)
{
   WriteVector(f,va,ldBinary);
   WriteFloat(f,&occ,1,ldBinary);
   if (!ldBinary) fprintf(f,"\n");
}

/* DumpTrAcc: dump transition acc to file f */
static void DumpTrAcc(FILE *f, TrAcc *ta, Boolean ldBinary)
{
   WriteMatrix(f,ta->tran,ldBinary);
   WriteVector(f,ta->occ,ldBinary);
   if (!ldBinary) fprintf(f,"\n");
}
FILE * baoDumpAccsParallel(HMMSet *hset, char *fname, int n, UPDSet uFlags, int index, DVector *noisyMeanData,DVector *noisyVarData,DVector *noisyMeanData2,DVector *noisyVarData2, Boolean ldBinary)
{
	FILE *f;
	HLink hmm;
	HMMScanState hss;
	int m,s;
	MixPDF* mp;

	fprintf(stderr, "\n# Output to %s, ldBinary = %d\n", fname, ldBinary);

	int vSize = hset->vecSize;

	Vector vTempNew = NULL;

	int numComponents = 0;
	vTempNew = CreateVector(&gstack,vSize);

	f = GetDumpFile(fname,n);

	NewHMMScan(hset, &hss);

	do {

		hmm = hss.hmm;
		DumpPName(f,hss.mac->id->name);
		WriteInt(f,(int *)&hmm->hook,1,ldBinary);

		while (GoNextState(&hss,TRUE)) {
			while (GoNextStream(&hss,TRUE)) {
				DumpWtAcc(f,((WtAcc *)hss.ste->hook)+index,outLdBinary);
				if (hss.isCont){
					while (GoNextMix(&hss,TRUE)) {
						numComponents++;
						MuAcc* tmpMu = NULL;
					  if ((uFlags&UPMEANS) && (!IsSeenV(hss.mp->mean))) {
							/* Thang v110701: mu */
						  tmpMu = ((MuAcc *)GetHook(hss.mp->mean))+index;
								int ii = 1;

								if(tmpMu->occ < 1e-40){ /* use non-compensated mu */
									DumpMuAccFromAscii(f,tmpMu->mu,tmpMu->occ,outLdBinary);
								} else { /* use compensated mu */
									for(; ii<=vSize; ii++){
										/*hss.mp->mean[ii] was updated when dump MMF, so It was processed using TMPC*/
										vTempNew[ii] = (noisyMeanData[numComponents][ii] - noisyMeanData2[numComponents][ii])*tmpMu->occ; /* thangLoad: vTemp[ii] = hss.mp->mean[ii] + tmpMu->mu[ii]/occ; */

										/* check if mu is infinite */
										if(!finite(vTempNew[ii])){
											fprintf(stderr, "# Mu infinite %d: %e\n", ii, vTempNew[ii]);
											exit(1);
										}

									}
									DumpMuAccFromAscii(f,vTempNew,tmpMu->occ,outLdBinary);
								}


						  TouchV(hss.mp->mean);
						}


						if ((uFlags&UPVARS) && (!IsSeenV(hss.mp->cov.var))) {
							/* Thang v110701: va */
							VaAcc * tmpVa = ((VaAcc *)GetHook(hss.mp->cov.var))+index;

								int ii = 1;
								if(tmpVa->occ < 1e-40 || tmpMu->occ < 1e-40){ /* use non-compensated var */
									if(tmpVa->occ > 0){
										fprintf(stderr, "#! Small muOcc=%e, varOcc=%e, use non-compensated var\n", tmpMu->occ, tmpVa->occ);
									}

									DumpVaAccFromAscii(f,tmpVa->cov.var,tmpVa->occ,outLdBinary);
								} else { /* Use compensated var */
									float muOcc = tmpMu->occ;
									for(; ii<=vSize; ii++){
										float ratio = vTempNew[ii]/muOcc; /* thangLoad: tmpMu->mu[ii]/muOcc; */
										vTempNew[ii] = (noisyVarData[numComponents][ii] + ratio*ratio)*tmpVa->occ;	/* thangLoad: vTemp[ii] = tmpVa->cov.var[ii]/occ - ratio*ratio */

										/* check if var is infinite */
										if(!finite(vTempNew[ii])){
											fprintf(stderr, "# Var infinite %d: %e\n", ii, vTempNew[ii]);
											exit(1);
										}
									}
									DumpVaAccFromAscii(f,vTempNew,tmpVa->occ,outLdBinary);
								}


							TouchV(hss.mp->cov.var);
						}
					}
				}
			}
		}
		if (!IsSeenV(hmm->transP)){
			DumpTrAcc(f, ((TrAcc *) GetHook(hmm->transP))+index,outLdBinary);
			TouchV(hmm->transP);
		}
		DumpMarker(f, outLdBinary);
	} while (GoNextHMM(&hss));
	EndHMMScan(&hss);



	/*fprintf(stderr, "# Done writing to %s\n", fname);*/

	return f;
}


int main(int argc, char *argv[])
{
	char*s;
	/* TPMC*/
	int numPoints, numObsVectors, isDebug;
	double alpha=0;
	char *option, *noiseAddOpt, *covarOpt, *wOpt;
	DMatrix dctDMatrix,cepDMatrix,invDctDMatrix, invCepDMatrix, Sq, W, VS;


	/* End TMPC*/

	void Initialise(char *hmmListFn);
	/*void UpdateModels(void);*/
	/*void StatReport(void);*/

	char newFn[MAXSTRLEN];
	FILE *f;
	Source src;
	float x; int i;
	char *filenameMMF;
	setbuf(stdout, NULL); /*unbuffered output.*/
	InitShell(argc,argv,dspr_version,dspr_vc_id);
	InitMem();
	InitMath();
	InitSigP();
	InitAudio();
	InitWave();
	InitVQ();
	InitModel();
	InitParm();
	InitLabel();
	InitTrain();
	InitUtil();
	InitFBLat();
	InitExactMPE();
	InitArc();
	InitDict();
	InitLat();
	InitNet();
	InitAdapt(&xfInfo);

	if (!InfoPrinted() && NumArgs() == 0)
		ReportUsage();
	if (NumArgs() == 0) Exit(EXITSTATUS);
	SetConfParms(); /*!Set up config variables.*/
	CreateHeap(&hmmStack,"HmmStore", MSTAK, 1, 1.0, 50000, 500000);
	CreateHMMSet(&hset,&hmmStack,TRUE);


	CreateHeap(&hmmStack2,"HmmStore", MSTAK, 1, 1.0, 50000, 500000);
	CreateHMMSet(&hset2,&hmmStack2,TRUE);

	while (NextArg() == SWITCHARG) {
		s = GetSwtArg();
		switch(s[0]){
			case 'a':
				if(!strcmp(s, "a")){
				  xfInfo.useInXForm = TRUE;
				} else if(!strcmp(s, "accs")){
					if (NextArg()!=STRINGARG)
						HError(-1,"dspr: accs prefix expected");
					accsPrefix = GetStrArg();
				} else if(!strcmp(s, "asciiin")){
					accsMode = ACCS_ASCII_IN;
					updateMode = UPMODE_DUMP; /* Dump acccs files with the updated means and variances specified by accsPrefix */
				} else if(!strcmp(s, "asciiout")){
					accsMode = ACCS_ASCII_OUT;
					updateMode = 0; /* no dump or update, simply output extraced means and variances in ASCII into accsPrefix files */
				}else HError(1, "Unknown option %s",s);

				break;

			case 'o':
							if (!strcmp(s, "o")) {
								if (NextArg()!=STRINGARG)
									HError(2319,"dspr: HMM file extension expected");
								newExt = GetStrArg();

							} else if (!strcmp(s, "outB")) {
								outLdBinary = TRUE;
							}else if(!strcmp(s, "option")) {
								option=GetStrArg();
							}
							break;

			case 'B':
							saveBinary=TRUE;
							break;
			case 'c':
							if(!strcmp(s,"covarOpt")){
								covarOpt = GetStrArg();
							}

							break;

			case 'H':
								if (NextArg() != STRINGARG)
									HError(2319,"dspr: HMM macro file name expected");
								if(!strcmp(s,"H")){
									char *x = GetStrArg();
									AddMMF(&hset,x);

								} else if(!strcmp(s,"Hprior")){
									if(!hset_prior_initialised){
										hset_prior_initialised=TRUE;
										CreateHMMSet(&hset_prior,&hmmStack,TRUE);
									}
									AddMMF(&hset_prior, GetStrArg());
								} else HError(1, "Unknown option %s", s);
								break;
			case 'i':
								if (!strcmp(s, "isDebug")) {
									isDebug =  GetChkedInt(0,1,s);
								}
								else if (!strcmp(s, "inB")) {
									inLdBinary = TRUE;
								}
								break;
			case 'n':
								if (!strcmp(s, "noiseAddOpt")) {
									noiseAddOpt = GetStrArg();
								}else if(!strcmp(s, "numPoints")){
									numPoints = GetChkedInt(0,100,s);
								}else if(!strcmp(s, "numObsVectors")){
									numObsVectors = GetChkedInt(0,100,s);
								} else HError(1, "Unknown option %s", s);

								break;

			case 'M':
								if (NextArg()!=STRINGARG)
									HError(2319,"dspr: Output macro file directory expected");
								newDir = GetStrArg();
								break;

			case 'P':
								if(!strcmp(s,"P")){
									char *x = GetStrArg();
									AddMMF(&hset2,x);

								}
								break;
			case 'T':

								trace = GetChkedInt(0,0100000,s);

								break;

			case 'w':
							if(!strcmp(s,"wOpt")){

								wOpt = GetStrArg();
							}

								break;


			default:
					HError(2319,"dspr: Unknown switch %s",s);
					break;
		}
	} /*matches while(NextArg() == SWITCHARG)*/

		InitVocab(&vocab); /* The actual dict is not needed, only the structure; this relates to HNet and reading lattices. */

	if (NextArg() != STRINGARG)
		HError(2319,"dspr: file name of hmm list expected");

	Initialise(GetStrArg()); /*GetStrArg() will return the hmmList.*/

	int numStates = hset.numStates;
	int numMix = hset.numMix;
	int vSize = hset.vecSize;

	DVector *stateCleanMeanData1 = malloc(numMix*sizeof(DVector));
	DVector *stateCleanVarData1 = malloc(numMix*sizeof(DVector));
	char **arCleanInfo1 = malloc((3*numStates+3)*sizeof(char*));
	DVector *stateCleanMeanData2 = malloc(numMix*sizeof(DVector));
	DVector *stateCleanVarData2 = malloc(numMix*sizeof(DVector));
	char **arCleanInfo2 = malloc((3*numStates+3)*sizeof(char*));
	DVector *stateCleanMeanData3 = malloc(numMix*sizeof(DVector));
	DVector *stateCleanVarData3 = malloc(numMix*sizeof(DVector));
	char **arCleanInfo3 = malloc((3*numStates+3)*sizeof(char*));

	DVector *stateCleanMeanData4 = malloc(numMix*sizeof(DVector));
	DVector *stateCleanVarData4 = malloc(numMix*sizeof(DVector));
	char **arCleanInfo4 = malloc((3*numStates+3)*sizeof(char*));


	loadMMF(&hset2,inLdBinary,stateCleanMeanData4,stateCleanVarData4,arCleanInfo4);   /*src is the still-open file.*/


	char accsPrefix1[256];
	char accsPrefix2[256];
	char accsPrefix3[256];
	char accsPrefix4[256];
	if(accsMode&ACCS_ASCII_OUT || accsMode&ACCS_ASCII_IN){
		sprintf(accsPrefix1, "%s%d", accsPrefix, 1);
		sprintf(accsPrefix2, "%s%d", accsPrefix, 2);
		sprintf(accsPrefix3, "%s%d", accsPrefix, 3);
		sprintf(accsPrefix4, "%s%d", accsPrefix, 4);
	}

	do {
		char *accfn;
		if(parMode == 0){  /* 0 means gather together the parallel files.*/

			/*check no of files:*/
			if(THREEACCS && MMIPrior){
				if((NumArgs() % 4) != 0)
					HError(1, "dspr: nAccs should divide into 4 for MPE with MMI Prior");
			} else if(THREEACCS){
				if((NumArgs() % 3) != 0)
					HError(1, "dspr: nAccs should divide into 3 for MEE&&MPEStoreML");
			} else if (!ML_MODE && (NumArgs() % 2) != 0)/*need an even number because .1 and .2 files.*/
				HError(2319,"dspr: Odd num of training files for single pass training");

			accfn = GetStrArg();

			/* Load first accumulator file (.1) */


			/* Warn if it is not a .acc.1 file */
			if(!strstr(accfn,".acc.1")) HError(-1, "Expecting a *.acc.1 file, got %s", accfn);

			src=thangLoadAccsParallel(&hset, accfn,  uFlagsAccs, 0, accsPrefix1,inLdBinary,stateCleanMeanData1,stateCleanVarData1,arCleanInfo1);   /*src is the still-open file.*/


			/*printf("10=%s",arCleanInfo1[10]);*/
			/*ShowDVector("",stateCleanVarData1[101],39);*/

			ReadFloat(&src,&x,1,TRUE);    /*x is the average log prob (MPE->avg correctness) */
			ReadInt(&src,&i,1,TRUE);      /* i the number of timeperiods (MPE->num correct words).*/
			CloseSource( &src );

			if(!MPE){ totalPr1 += x*i;  totalT += i; }
			else {
				TotalCorr += x*i; TotalNWords += i;
			}
			/*fprintf(stderr, "# TotalCorr=%e, TotalNWords=%d, x=%e, i=%d\n", TotalCorr, TotalNWords, x, i);*/

			/* Load second accumulator file (.2) */
			if(!ML_MODE){                   /* ..Then load MMI acc. */
				accfn = GetStrArg();

				if(!strstr(accfn,".acc.2")) HError(-1, "Expecting a *.acc.2 file, got %s", accfn);

				src=thangLoadAccsParallel(&hset, accfn,  uFlagsAccs, 1, accsPrefix2,inLdBinary,stateCleanMeanData2,stateCleanVarData2,arCleanInfo2);

				/*printf("arCleanInfo1[101]=%s",arCleanInfo2[101]);*/
				ReadFloat(&src,&x,1,TRUE); /*!x must be the average log prob*/
				ReadInt(&src,&i,1,TRUE); /*!and i the number of timeperiods.*/
				CloseSource( &src );
				/*pr2 contains the MMI stats.*/
				totalPr2 += x*i;  /*totalT += i;*/
				/*fprintf(stderr, "# totalPr2=%e, x=%e, i=%d\n", totalPr2, x, i);*/
			}

			/* Load third accumulator file (.3) */

			if(THREEACCS){             /* (MPE case). */
				accfn = GetStrArg();

				if(!strstr(accfn, ".3")) HError(1, "Error, expecting a HDR?.acc.3 file, got %s", accfn);
				src=thangLoadAccsParallel(&hset, accfn, uFlagsAccs, 2/*third position(1st==0) of hset*/, accsPrefix3,inLdBinary,stateCleanMeanData3,stateCleanVarData3,arCleanInfo3); /*src is the still-open file.*/

				ReadFloat(&src,&x,1,TRUE); /*!x must be the average log prob*/
				ReadInt(&src,&i,1,TRUE); /*!and i the number of timeperiods.*/
				CloseSource( &src );
				totalPr1 += x*i;  totalT += i;
			}


		} /*[parMode]*/
	} while (NumArgs()>0);

updateMode = UPMODE_DUMP;
/* Start Saving the statistics to  HDR format*/
	if (updateMode&UPMODE_DUMP ){

		MakeFN("noisyHDR$.acc.1",newDir,NULL,newFn);

		f=baoDumpAccsParallel(&hset,newFn,parMode, uFlagsAccs, 0, stateCleanMeanData1,stateCleanVarData1, stateCleanMeanData4,stateCleanVarData4,outLdBinary);
		if(MPE){ /* .1 acc contains the MPE crit at the end. */
			x = TotalCorr/TotalNWords;
			WriteFloat(f, &x, 1, outLdBinary);
			WriteInt(f, &TotalNWords, 1, outLdBinary);
		}  else {
			x = totalPr1/totalT;
			WriteFloat(f,&x,1,outLdBinary);
			WriteInt(f,&totalT,1,outLdBinary);
		}
		fclose( f );
		puts("Done Saving to noisyHDR0.acc.1");
		if(!ML_MODE){ /* the dump .2 accs. */
			MakeFN("noisyHDR$.acc.2",newDir,NULL,newFn);
			f=baoDumpAccsParallel(&hset,newFn,parMode, uFlagsAccs, 1, stateCleanMeanData2,stateCleanVarData2, stateCleanMeanData4,stateCleanVarData4,outLdBinary);
			x = totalPr2/totalT; /*MMI den prob in either MMI or MPE case.*/
			WriteFloat(f,&x,1,outLdBinary);
			WriteInt(f,&totalT,1,outLdBinary);
			fclose( f );

			puts("Done Saving to noisyHDR0.acc.2");
		}

		if(THREEACCS){   /* third acc for mle smoothing of MPE. */
			MakeFN("noisyHDR$.acc.3",newDir,NULL,newFn);
			f=baoDumpAccsParallel(&hset,newFn,parMode, uFlagsAccs, 2, stateCleanMeanData3,stateCleanVarData3, stateCleanMeanData4,stateCleanVarData4,outLdBinary);
			x = totalPr1/totalT; /* This is where the MLE prob is stored. */
			WriteFloat(f, &x, 1, outLdBinary);
			WriteInt(f, &totalT, 1, outLdBinary);
			fclose( f );
			puts("Done Saving to noisyHDR0.acc.3");
		}

		if(trace&T_TOP) PrintCriteria();
	}
	if (parMode <= 0){      /*parMode <= 0, so do re-estimation.*/
		if (stats) StatReport();
		if (uFlags&UPXFORM) {/* ensure final speaker correctly handled */
			UpdateSpkrStats(&hset,&xfInfo, NULL);
			if (trace&T_TOP) {
				printf("Reestimation complete - average log prob per frame = %e (%d frames)\n",
						totalPr1/totalT, totalT);
			}
		}

	}
	Exit(EXITSTATUS);
	return (0);          /* keep compiler happy */
}

/* -------------------------- Initialisation ----------------------- */

void Initialise(char *hmmListFn)
{
	char buf[256];

	   CreateHeap(&transStack,   "transStore",    MSTAK, 1, 0.5, 1000,  10000);
	   CreateHeap(&accStack,   "accStore",    MSTAK, 1, 1.0, 50000,  500000);
	   CreateHeap(&latStack,"latStore", MSTAK, 1, 1.0, 50000, 500000);


	   /* Load HMMs and init HMMSet related global variables */


	   MakeHMMSet( &hset, hmmListFn );
	   LoadHMMSet( &hset,hmmDir,hmmExt);
	   SetParmHMMSet(&hset);

	   MakeHMMSet( &hset2, hmmListFn );
	   LoadHMMSet( &hset2,hmmDir,hmmExt);
	   SetParmHMMSet(&hset2);



	   if(hset_prior_initialised){
	     MakeHMMSet( &hset_prior, hmmListFn );
	     LoadHMMSet( &hset_prior, hmmDir,hmmExt);
	   }

	   if(MMIPrior) NumAccs=4;
	   else if(ML_MODE) NumAccs=1;
	   else if(THREEACCS/*MPE||MPEStoreML*/) NumAccs=3;
	   else NumAccs=2;

	   {
	      uFlagsAccs =  uFlags|(uFlags&UPMEANS||uFlags&UPVARS ? UPMEANS|UPVARS : 0);
	      /*That modification to uFlags means: if either mean or var is updated, accumulate both.*/
	      AttachAccsParallel(&hset, &accStack, uFlagsAccs, NumAccs);
	      ZeroAccsParallel(&hset, uFlagsAccs, NumAccs);
	   }


	   P = hset.numPhyHMM;
	   L = hset.numLogHMM;
	   vSize = hset.vecSize;
	   S = hset.swidth[0];
	   hsKind = hset.hsKind;

	   if(S > 1)
	      HError(-1, "dspr: Code is intended to support multiple streams but code has not been debugged.  Be warned!");
	   /*Check that hset is the right kind.*/
	   if(!(hsKind == PLAINHS || hsKind == SHAREDHS || hsKind == TIEDHS))
	      HError(1, "dspr: hset kind not PLAIN or SHARED or TIED.");
	   if((!(hset.ckind == NULLC)) && (!(hset.ckind == DIAGC)))
	      HError(1, "dspr: cov kind not DIAGC.");

	   /* Additional code for the adaptation updates */
	   /*!Deleted*/

	   /*Initialise those modules.*/
	   InitialiseFBInfo(&fbInfo, &hset,   uFlags|(uFlags&UPMEANS||uFlags&UPVARS ? UPMEANS|UPVARS : 0), twoDataFiles);
	   /*That modification to uFlags means: if either mean or var is updated, accumulate both.*/


	   /* Set the variance floor */
	   SetVFloor( &hset, vFloor, minVar); /*This sets the array minVar up... but it only works if
	                                        minimum variances have been set within the hmmset (e.g,
	                                        by HCompV.*/

	   if (trace&T_TOP) {
	      printf("HAccsConvert  Updating: ");
	      if (uFlags&UPTRANS) printf("Transitions ");
	      if (uFlags&UPMEANS) printf("Means ");
	      if (uFlags&UPVARS)  printf("Variances ");
	      if (uFlags&UPMIXES && maxM>1)  printf("MixWeights ");
	      if (uFlags&UPXFORM)  {
	        if (xfInfo.inSpkrPat == NULL) xfInfo.inSpkrPat = xfInfo.outSpkrPat;
	        if (xfInfo.paSpkrPat == NULL) xfInfo.paSpkrPat = xfInfo.outSpkrPat;
	        if (uFlags != UPXFORM)
	          HError(999,"Can only update linear transforms OR model parameters!");
	        xfInfo.useOutXForm = TRUE;
	        /* This initialises things - temporary hack - THINK!! */
	        CreateAdaptXForm(&hset, "tmp");
	      }
	      printf("\n ");
	      if (parMode>=0) printf("Parallel-Mode[%d] ",parMode);
	      printf(" - covkind is %s; ",CovKind2Str(hset.ckind,buf));
	      printf(" system is ");
	      switch (hsKind){
	      case PLAINHS:  printf("PLAIN\n");  break;
	      case SHAREDHS: printf("SHARED\n"); break;
	      case TIEDHS:   printf("TIED\n");   break;
	      case DISCRETEHS: printf("DISCRETE\n"); break;
	      }

	    /*  printf("%d Logical/%d Physical Models Loaded, VecSize=%d\n",L,P,vSize);*/
	      if (hset.numFiles>0)
	         printf("%d MMF input files\n",hset.numFiles);
	      if (mmfFn != NULL)
	         printf("Output to MMF file:  %s\n",mmfFn);
	      fflush(stdout);
	   }

	   if (parMode != 0){
	      ConvDiagC(&hset,TRUE); /*!Converts DIAGC to INVDIAGC: changing over the CovKind in each HMM and inverting all the elements.*/
	      ConvLogWt(&hset);
	   }
}

/* ------------------- Statistics Reporting  -------------------- */

/* PrintStats: for given hmm */
void PrintStats(FILE *f, int n, HLink hmm, int numEgs)
{
   WtAcc *wa;
   char buf[MAXSTRLEN];
   StateInfo *si;
   int i,N;

   N = hmm->numStates;
   ReWriteString(HMMPhysName(&hset,hmm),buf,DBL_QUOTE);
   fprintf(f,"%4d %14s %4d ",n,buf,numEgs);
   for (i=2;i<N;i++) {
      si = hmm->svec[i].info;
      wa = (WtAcc *)((si->pdf+1)->hook);
      fprintf(f," %10f",wa->occ);
   }
   fprintf(f,"\n");
}

/* StatReport: print statistics report */
void StatReport(void) /*This is used by other programs so I have had to change it back to how it was.*/
{
   HMMScanState hss;
   HLink hmm;
   FILE *f;
   int px;

   if ((f = fopen(statFN,"w")) == NULL){
      HError(2311,"StatReport: Unable to open stats file %s",statFN);
      return;
   }
   NewHMMScan(&hset,&hss);
   px=1;
   do {
      hmm = hss.hmm;
      PrintStats(f,px,hmm,(int)hmm->hook);
      px++;
   } while (GoNextHMM(&hss));
   EndHMMScan(&hss);
   fclose(f);
}


/* UpdateModels: update all models and save them in newDir if set,
   new files have newExt if set */
void baoUpdateModels(DVector *mean, DVector *var, int numFeatures)
{
   HMMScanState hss;

   int i;
   {     /*Measure total occupancy and update Gaussians. */
      int k=0;
      if(hset.hsKind==PLAINHS || hset.hsKind==SHAREDHS){
         NewHMMScan(&hset, &hss);
         while(GoNextMix(&hss, FALSE)){ /* Update mixes. */
            k++;

            for(i=1;i<=numFeatures;i++)
            {
            	hss.mp->mean[i] = mean[k][i];
            	hss.mp->cov.var[i] = var[k][i];
            }
            /*FixGConsts(hss.hmm);*/
            FixDiagGConst(hss.mp);
            /*hss.mp->cov.var = (SVector)var[k];*/

           /*if(UpdateGauss(hss.s, hss.mp)) nMix++; else nLeft++;*/
         }
         EndHMMScan(&hss);
      }
   }


     if (trace&T_TOP) {

      if (mmfFn == NULL)
         printf("Saving hmm's to dir %s\n",(newDir==NULL)?"Current":newDir);
      else
         printf("Saving hmm's to MMF %s\n",mmfFn);
      fflush(stdout);
   }
   ClearSeenFlags(&hset,CLR_ALL); /*?*/

   SaveHMMSet(&hset,newDir,newExt,NULL,saveBinary);
   if (trace&T_TOP) {
      printf("Saved hmm\'s\n");
      PrintCriteria();
   }
}


/* ----------------------------------------------------------- */
/*                      END:  dspr.c                       */
/* ----------------------------------------------------------- */
